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Abstract. We reexamine the interpretation of the annual modulation signal observed by 
the DAMA experiment as due to WIMPs with a spin-dependent coupling mostly to protons. 
We consider both axial-vector and pseudo-scalar couplings, and elastic as well as endothermic 
and exothermic inelastic scattering. We conclude that the DAMA signal is in strong tension 
with null results of other direct detection experiments, particularly PICASSO and KIMS. 
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1 Introduction 


Direct Dark Matter (DM) searches aim at detecting Weakly Interacting Massive Particles 
(WIMPs) scattering with nuclei in a target material. The most stringent limits on the DM 
mass and cross section parameter space are currently set by LUX ppj and SuperCDMS [2] for 
WIMPs with spin-independent interactions and spin-dependent interactions with neutrons, 
and by PICASSO [3], SIMPLE [3], COUPP [5], and KIMS [B] for those with protons. While 
none of these experiments have detected a DM signal so far, three other experiments have 
signals that can be interpreted as due to WIMP scattering: DAMA and DAMA/LIBRA 
(called here DAMA from now on) (7) , CoGeNT [HI E] , and CDMS-II with silicon detectors 
[10]. CRESST-II [TT] has not confirmed a previous DM hint found by the same collaboration 
[12]. Among the potential DM signals, DAMA’s observation of an annually modulated rate 
has the highest statistical significance. However, these potential signals are challenged by 
the null results of other direct detection experiments which exclude the possibility of WIMP 
scattering in a large number of particle models. In particular, the scattering cross section 
fitting the DAMA data for WIMPs with isospin-conserving spin-independent interactions is 
several orders of magnitude above the 90% CL LUX limit SSI- 

Over the years, several particle candidates have been proposed with properties which 
enhance a potential signal in DAMA and weaken the main limits imposed by direct DM 
searches with negative results, among them WIMPs with spin-dependent coupling mostly to 
protons, which we reevaluate in this paper. 

Inelastic DM [141420] scatters to another particle state, either heavier (endothermic 
scattering) or lighter (exothermic scattering, see e.g. [211 l22j ). when colliding with a nucleus. 
Endothermic scattering favors heavier targets, thus enhancing scattering off I in DAMA while 
reducing scattering off lighter targets such as Ge. Moreover, this type of interaction enhances 
the annual modulation amplitude, thus pushing the cross section needed to fit the DAMA 
data to lower values. However, experiments employing Xe as target material, which is heavier 
than I, rule out endothermic scattering of DM as an explanation to the DAMA data unless 
there is an additional feature of the interaction that favors a signal in DAMA. Two types 
of WIMP couplings favor Na and I (DAMA) over Xe (LUX, XENONIO) and Ge (CDMS, 
SuperCDMS): a spin-dependent coupling mostly to protons and a magnetic dipole moment 
coupling [191123 ]. The reason for the first is that the spin of a nucleus is mostly due to an 
unpaired nucleon and Na and I have an unpaired proton, while Xe and Ge have an unpaired 
neutron. The reason for the second is the large magnetic moment of both Na and I. We 
do not consider here inelastic magnetic DM, which Ref. |20| recently found still marginally 
compatible with all negative results of direct searches, for a small value of the I quenching 
factor (without the aid of inelasticity, instead, this candidate has been shown to be ruled out 

EH). 

The possibility that a DM candidate with spin-dependent interactions mostly with pro¬ 
tons would explain the DAMA signal was, to the best of our knowledge, first studied in 
Ref. [25] . Compared to spin-independent interactions, spin-dependent couplings reduce the 
bounds from experiments with heavy targets, most notably LUX, due to the lack of the 
usual Ay enhancement factor proper to the spin-independent interaction (At being the mass 
number of the target nucleus). This interaction might explain why the DAMA signal is not 
seen by LUX and SuperCDMS. In this case, bounds from PICASSO, SIMPLE, COUPP, and 
KIMS become relevant, since they contain F, I and Cs, all nuclei with an unpaired proton. 
This candidate was further studied in Ref. [26], in the context of both elastic and inelastic 
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endothermic scatterings. The inelastic endothermic kinematics reduces the expected rate in 
experiments employing F (PICASSO, SIMPLE) because it is light, thus making the COUPP 
(CF 3 I) and KIMS (Csl) bounds the most relevant constraints on WIMP scatterings off I 
in DAM A. Ref. |2Dj found that a small portion of the parameter space favored by DAM A 
for inelastic spin-dependent couplings with protons can still escape all bounds from null 
experiments. 

Inelastic exothermic scattering EH [27] favors lighter targets, so it favors Na in DAM A 
over heavier nuclei (Ge and Xe). In this case the most important limits are set by experiments 
containing F (PICASSO and SIMPLE). 

Recently, Ref. [28] studied a Dirac WIMP candidate coupled to standard model (SM) 
fermions through a light pseudo-scalar mediator, and claimed that with a contact interaction 
and elastic scattering it reconciles the DAMA data with the null results of other experiments 
at the 99% credible level. The model produces a non-standard spin-dependent interaction, 
with the noteworthy feature that, for universal flavor-diagonal quark couplings to the pseudo¬ 
scalar mediator, the WIMP couples mainly to protons. The couplings of pseudo-scalar light 
bosons (rritf, < 7 GeV) with quarks are strongly constrained by rare meson decays [291431] . 
and unless the pseudo-scalar coupling to the DM (called #dm below) is very large, <?dm > 10 3 , 
the one-particle exchange scattering cross section required in Ref. [28] is rejected [31). The 
flavor physics bounds on pseudo-scalar couplings to quarks proportional to the quark mass 
are less stringent EH, but in this case Ref. [28] found that the resulting proton to neutron 
coupling ratio is not large enough to reconcile a DM signal in DAMA with the results of 
other direct detection experiments. Leaving aside the limits from other types of experiments 
we concentrate here on direct detection. 

Ref. [28] employed a Bayesian analysis, where a number of uncertain parameters like 
quenching factors and background levels, as well as astrophysical quantities, are marginalized 
over. While the process of marginalization (i.e. integrating over nuisance parameters with 
assumed prior probability distributions) is the proper treatment of uncertain and uninterest¬ 
ing parameters in the context of Bayesian statistics, it makes it unclear whether there exists 
at least one set of values of the uncertain parameters, in particular one halo model, that 
produces the same result of the analysis. 

In this paper we reconsider the viability of a signal due to WIMPs with spin-dependent 
coupling mostly to protons as an explanation of the DAMA data. We study both axial-vector 
and pseudo-scalar couplings, which lead respectively to s x ■ s p and (s x -q){s p -q) couplings in 
the non-relativistic limit (s x and s p are the spins of the WIMP x and the proton respectively, 
and q is the momentum transfer). We assume the mediator to be either heavy enough for the 
contact interaction limit to be valid, or otherwise much lighter than the typical momentum 
transfer of the scattering process than the typical momentum transfer of the scattering process 
(we refer to this last case as “massless”). The possibilities of elastic and inelastic scattering, 
both endothermic and exothermic, are considered. 

In Section [2] we present the differential cross sections for axial-vector and pseudo-scalar 
couplings, which can be used in the direct detection rate formula in Section [3] The analysis 
methods we adopt for experimental data are described in Section [I] and our results assuming 
a standard model of the dark halo of our galaxy are presented in Section [5] In Section [ 6 ] we 
describe our halo-independent analysis and present the related results. Our conclusions are 
given in Section [7] 
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2 Cross sections 


2.1 Axial-vector (AV) interaction 

An AV coupling leads to the usual spin-dependent interaction. The effective Lagrangian for 
the elastic scattering of a DM particle x, either a Dirac or a Majorana fermion, with AV 
couplings to nucleons, mediated by a vector boson of mass m^, is 


« 2 av 


ffPM 

2(mj ~ q»q P ) 


22 . 

N=p,n 


( 2 . 1 ) 


Here we assumed a one-particle exchange process, q^ is the momentum transfer four-vector, 
and A is a nucleon, p or n. 5 dm and ojv are the mediator coupling constants to x and N, 
respectively, and they are real. The scattering amplitude is 


./#AV 


gPM 

2 (ml - q^q^) 


a N u x Yl U X U N ~ 1 ^ u N . 

N=p,n 


( 2 . 2 ) 


We now follow Ref. [32] because we will largely use the nuclear form factors given in this 
reference. We first take the non-relativistic limit of the Dirac spinors, in the chiral represen¬ 
tation, for both x and N: u s (p) ~ ^/l/4m((2m — p ■ 8) £ s , (2m + p ■ 8) £ s ) , where 8 are the 
Pauli matrices. This limit is justified by the fact that the DM initial speed and the exchanged 
momentum are small. The matrix element for scattering off a single nucleon then reads 


= —8 m^m 


5pm 

ml + q 2 


22 aN ■ (+v) > 

N=p,n 


(2.3) 


where (s x ) = an d (+v) = ( see e -§- Eqs. (44), (47d), and (49) of Ref. (33]). 

Notice that this matrix element assumes the usual form for the Dirac spinors with normaliza¬ 
tion u s (p)u s {p) = 2 m5 ss ; Quantum Mechanical amplitudes usually assume a different state 
normalization, which differs by a factor of 2y/m 2 + p' 2 . With this normalization Eq. (2.3) 
would be replaced by ,£^2 = “ 2 5pm (q 2 + ^) _1 Y,N= P) n a N (s x ) ' (sn)- 

For a model of inelastic DM, one could introduce two Dirac fields, xi and X2, with 
slightly different masses and the DM-nucleon effective Lagrangian 


-£av = 


5pm 


2 (ml + q 2 ) 


22 aN X27 M 7 5 Xi n X p 1 5 N + h.c. 


(2.4) 


N=p,n 


The ^dm coupling can now be complex. We use the same symbol, 5 dm, for the couplings in 


Eqs. (|2.1|) and (|2.4|), because then the expression of in Eq. (|2.6|) is valid both for elastic 


and inelastic scattering, xi is the DM particle entering the scattering process, with mass rn, 
while X '2 is the DM particle in the final state, with mass m! = m + 6. The sign of the mass 
splitting 5 determines the different kinematic regimes: 5 > 0 implies endothermic scattering, 
(5 < 0 implies exothermic scattering, while <5 = 0 implies elastic scattering. 

One may attempt to build an inelastic DM model without introducing additional degrees 


of freedom by assuming the interaction in Eq. (2.1) and adding a small Majorana mass 


term which produces two almost degenerate Majorana fermions, xi and X' 2 , in which case 
X = Xi +7X2 becomes a quasi-Dirac fermion. However, as noted in Ref. [26], this interaction 
produces diagonal terms which result in elastic scattering rather than inelastic. In this case, 
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one can instead write an effective tensor interaction which produces inelastic 

scattering since the diagonal interaction terms vanish identically. The non-relativistic limit 
of this operator is also s x - sn- 

The differential cross section for DM-nucleus scattering, for both the elastic and inelastic 


interactions introduced above in Eqs. (2.1) and (2.4) 


is 


dcr^ v 

d E^ 


= a. 


AV r ^T_ 

2 4 


m: 


+ 2 tutEr I v 


~o F av(q 2 ) 


(2.5) 


where Er = q 2 /2rriT is the nuclear recoil energy, v is the incoming WIMP speed, and g, p is 
the DM-proton reduced mass. F 2 y (q 2 ) is a nuclear form factor including spin dependence, 
and will be defined in Section 


2.3 


cr,^ v is the total DM-proton cross section in the limit of 


contact interaction 3> q = ^j2rrvpE^ 


a™ = 


3|5DM| 2 Op /x; 


47r 


m 


( 2 . 6 ) 


The term in parenthesis in Eq. (2.5) accounts for long-range interactions, when m$ is smaller 
than or comparable to q = \Z2 i71tEr. For typical target masses of a few tens of GeV and 
recoil energies around few to tens of keV, the interaction becomes effectively long-range if the 
mediator mass is smaller than several MeV: m<p <C q — 20 MeV^(E r/IO keV)(?nR/20 GeV). 

In order to plot our results for long-range interactions, we express the differential cross 
section in terms of a reference total cross section cr^ V,lcf = ffp V (m^, = m^ ef ) corresponding 
to a reference mediator mass ml, which we set equal to 1 GeV: 


dtJ T V = a AV, ref ™T ( 

dE R p 2g 2 + 2m T E R 


^2 F AviQ 2 ) ■ 


(2.7) 


The massless mediator limit thus corresponds to setting m^ = 0 in the equation above. In the 
following we will refer to any scenario with m ^ -C q 2 as massless mediator limit or long-range 
limit. 

Given that a large value of a p /a n is needed to suppress the strong LUX and SuperCDMS 
constraints, we will assume the maximally isospin-violating coupling a n = 0 for the AV 
interaction. 


2.2 Pseudo-scalar (PS) interaction 

Here the DM particle is a Dirac fermion y, coupled to a real PS boson 4> with mass m^, 

=2dm = ( 2 -8) 

(as in Refs. [28, [33]), with a real coupling constant 5 dm- The PS field couples also to the SM 
quarks with real coupling g q , 

■Z’q = ^gq&q^q ■ ( 2 -9) 

^ q 

While PS couplings to quarks are usually taken to be proportional to the fermion mass (see 
e.g. [33]), we will assume instead a flavor-universal coupling g q = g, which introduces a larger 
\a p /a n \ ratio, a p /a n — —16.4 [28] (see below). 
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To model inelastic scattering we assume a non-diagonal coupling of two Dirac DM fields 
Xi and X 2 with <j>, 

( 2 - 10 ) 


(fz? • ^DM / — 5 - -l 

^DM = -l—f=-<px 27 Xl + h.C. . 


V 2 

Again, the +dm coupling can now be complex. With this definition of c/dm , Eqs. ( |2.8[ ) and 
(2.10) yield the same expression for in Eq. (2.16). Diagonal interaction terms as well as 


non-diagonal mass terms can be forbidden by assuming a Z 2 symmetry under which both (j), 
Xi ( or X 2 ), and all the SM electroweak doublets have charge —1. 

The DM-nucleon effective Lagrangian, for elastic scattering (for inelastic scattering we 
should have X27 a 1 instead of X7 5 x)> and assuming one-particle exchange, is 


^PS = 


9dm 




2(771^ — q^Qn) J—* 
v <P y N=p,i 


aN 5 N 


and yields in the non-relativistic limit 


^PS = 2 


9dm 
m 2 , + q 2 


°ov((4> ■Q)i(sN)-q) , 


( 2 . 11 ) 


( 2 . 12 ) 


N=p,n 


where (s x ) and (sjv) were defined after Eq. (2.3). This is a different type of spin-dependent 


interaction than in Eq. (2.3). Due to the extra factors of q, the PS cross section receives a large 
q^/mj^m 2 suppression with respect to the AV cross section. Therefore, the normalization 
of the signal and its spectrum, and also the nuclear form factors are different in the two 
cases 


(see Section 2.3). Given the large momentum suppression, one needs to check the 
existence of unsuppressed radiative corrections to this tree-level cross section, that would 


spoil the setup. The PS interaction in Eq. (2.11) with a Dirac fermion x has been proven 
not to produce such corrections 
fermion. 


j, while this would not be the case if x were a Majorana 
The proton and neutron couplings appearing in Eq. (|2.11[) are given by 


a N = 9 ^2 

f=u,d,s 


m N 
m f l 


1 - 


E 


777. 


m r 


a r-> 


E 

f=u,d,s 


m N 
nrif 


m 


m f J 


AN) 


(2.13) 


with m = (1 jm u + 1/m^ + 1 /m s ) 1 . The subscripts / in (2.13) indicate quark flavors. The 
Ay) V ' ) factors parametrize the quark spin content of the nucleon, and are usually determined 
experimentally or computed with lattice calculations. As in Ref. [28j, we adopt the following 
values from 


a (p) = A^ n) = +0.84 , = Al n) = -0.44 , A<f } = A( n) = -0.03 , (2.14) 

with which a p — —OAg. As a natural feature of this model, the proton coupling a p is larger 
(in modulus) than the neutron coupling a n , by an amount that depends on the choice of 


, the values in Eq. (2.14) are conservative in the sense that 


the Af ! s. As noted in 
they minimize the ratio a p /a n with respect to other values encountered in the literature (see 
e.g. Table 4 in [33]). In this case a p /a n = —16.4. 

The differential cross section for the PS interaction, for both elastic and inelastic scat¬ 
tering, is 


d(jy S PS 3m|P^ 


1 




1 


d E r p 8+6 


,ref 4 l J7i2 _|_ 2 TIItEr I V 


4e| s (9 2 ) • 


(2.15) 
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Fp S (q 2 ), to be defined in Section 2.3, is the nuclear form factor including spin dependence. 


PS\ 

<7p S is the total DM-proton cross section in the limit of contact interaction, 


^ S = 


IsdmI /4 u ref4 


12tt 


m ^ m 2 m 2 


(2.16) 


In this case, the total DM-proton cross section has a u 4 dependence, and, for the purpose of 
plotting our results in terms of the reference cross section, in Eq. (2.16) we evaluate cr 4>s 


at a 

reference speed v rel . We set u ref equal to the rotational speed of our Local Standard of Rest, 
220 krn/s, which is representative of the WIMP speeds with respect to Earth. 

For long-range PS interactions we proceed in the same manner as for the long-range 
AV interactions, by writing the differential cross section in terms of a reference total cross 

PS, ref _ PS 


section <7p ,re = (m^ = m^ ei ), with m™* = 1 GeV: 


du? 3 

d£ R 


= a. 


PS, ref 31 




8 /4 


, 1 

-yref 4 \m^ + 2m T E R J v 


s(r) 


(2.17) 


2.3 Nuclear form factors 

We adopt the form factors computed in Ref. [32] using standard shell model techniques, for 
the nuclides for which they are available, namely the main stable isotopes of Ge, Xe, Na, I, 
and F. In these cases, we define 

F Av(Q 2 ) = ^2 E ^NaN'(F^’ N '\q 2 ) + F^ N '\q 2 )) (2.18) 

P N,N'=p,n 


for the AV interaction, and 

Fpsio^) = ~2 E a N aN'F^ N,) (q 2 ) (2.19) 

P N,N'=p,n 


for the PS interaction. The (squared) nuclear form factors EVy and EVy/ are tabulated in 
[32] for the nuclides mentioned above. These form factors can be employed unmodified also 
for inelastic scattering [20]. We include a factor of 1/3 in the definition of E/ v in order to 
normalize F \ v and Fp S to be 1 in the limit of zero momentum transfer when the target is an 
isolated proton. This factor traces back to EVv being twice as large as Es" at q = 0, which 
is consistent with the fact that EV// corresponds to the component of the nucleon spin along 
the direction of the momentum transfer, while Eyy corresponds to the transverse component. 

E^ v (g 2 ) can be expressed (see Eqs. (59), (60) and (77c) in Ref. [32j) in terms of the 
usual nuclear spin structure function S(q 2 ) = aQSoo(q 2 ) + aoaiSoi(q 2 ) + a 2 Sn(q 2 ) [37] (with 
o-o = cl p + a n and a\ = a p — a n the isoscalar and isovector parameters): 


F 2 av(q 2 ) = 


4-7T 


3(2 J T + 1) al 


S(q 2 


( 2 . 20 ) 


with Jt the spin of the target nucleus. At zero momentum transfer 


S( 0) 


1 (2 Jt + 1 )(Jt + 1) 

7r Jt 


(a p (S p ) + a n (S n )) 2 , 


( 2 . 21 ) 
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where (S p ) = = Jt\S p \Jt, Mt = Jt), and where S p is the component of S p = 

T protons along the z-axis [38j ((S n ) is defined analogously). Notice that (S p ) and (S n ) 
are often denoted with boldface style in the literature, although they are not vector quan¬ 
tities. v can then be expressed in terms of the usually called spin-dependent form factor 
^sW) = S(q 2 )/S(0) as 

f L(Q 2) = 4( ^ 1} (V P > + ^(Sn)y Fl D (q 2 ) . ( 2 . 22 ) 


For the nuclides for which no form factors have been computed in Ref. 


(Cl, C 


and Cs), we define F‘^ v (q 2 ) by means of Eq. (2.22), with the spin-dependent form factor 


m 


Gaussian form 


PsW) = S2/4 ; 


(2.23) 


here we take R = ^0.92A^/ 3 + 2.68 — 0.78^(A^/ 3 — 3.8) 2 + 0.2^ fm, with At the mass num- 

ber of the target nucleus [39]. In this case we also assume Ep S = Flv, which we expect to 
be approximately valid at low q 2 . For Cs, a component of KIMS’s target material, we take 
(S p ) = -0.370, (S n ) = 0.003 MM- For SIMPLE, we use (S p ) = -0.051, (S n ) = -0.0088 
for both 35 C1 and 37 C1 [38], and (S p ) = —0.026, (S n ) = —0.155 for 13 C j32j. Notice that there 
are large uncertainties in the hadronic matrix elements (S p ) and (S n ) and the nuclear form 
factors, which differ in different nuclear models (see e.g. Fig. 1 of Ref. |43| for the Xe nuclear 
structure functions). Factors of 2 difference in different calculations are not uncommon [44] . 


3 Direct detection rate 


The DM-nucleus scattering rate for a target nuclide T with mass tut is 


d Rt 
d Er ' 


(Er, t) = — f d 3 vvf(v,t)^F(E R ,v) 

m Jv>v min (E R ) dE R 


(3.1) 


where m is the WIMP mass, p is the local DM density, f(v, t ) is the DM velocity distribu¬ 
tion in Earth’s reference frame, and dax/dER is the differential cross section for a WIMP 
scattering off the target T. v m ^ (E- r) is the minimum WIMP speed needed to impart to the 
target nucleus a recoil energy Er. 

The detectors do not measure directly the recoil energy, but they measure a related 
energy E' (sometimes expressed in ‘keV electron equivalent’ or keVee, or else in number of 
photoelectrons). The rate in Eq. (3.1) is related to the event rate measured by an experiment 


within a detected energy interval [E \, E! 2 ] as 


R 




, w = E^f 

“ m T Jo 




dE' e(E', Er)G t (Er, E') 


(3.2) 


where Ct is the mass fraction of the nuclide T in the detector, Gt(Er, E') is the target- 
dependent resolution function of the detector and e(E', Er) is the experimental acceptance. 
The resolution function is defined as the probability distribution for a nuclear recoil energy 
Er to be measured as E', and incorporates the mean value (E') = Qt(Er)Er, with Qt the 
target’s quenching factor, and the detector energy resolution. The detector acceptance is a 
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function of both Er and E' in general, but it is often given to be Er independent, or the 
dependencies on Er and E' are factorizable. 

Except for our halo-independent analysis, in this paper we assume the Standard Halo 
Model (SHM) for the dark halo of our galaxy, where the DM local density is p = 0.3 
GeV/cm 3 and the velocity distribution of WIMPs in the galactic frame is a truncated 
Maxwell-Boltzmann distribution 

fG(u) = 1 exp(-U 2 /^o) %esc - u). (3.3) 

N esc [voV ^) 


Here vq is the velocity dispersion and v esc is the escape speed from our galaxy. The normal¬ 
ization factor 


AW = 


= erf(u esc /uo) - 2(v esc /vo) exp(—v 2 sc /uo)/\/vr (3.4) 

ensures that f d 3 uf G (u ) = 1. We consider vq to be the same as the velocity of the Local 
Standard of Rest vq = 220 km/s, and we take v esc = 533 km/s, according to recent Radial 
Velocity Experiment (RAVE) 2013 results |45j. The velocity distribution in Earth’s frame, 


f(v,t ) in Eq. (3.1), can be obtained with the Galilean transformation 


f(v , t) = f G (v & + v e (t) + v) , 


(3.5) 


where u© and u©(t) are the velocity of the Sun with respect to the galaxy and the time 
dependent velocity of Earth with respect to the Sun, respectively. We take w© = 232 km/s, 
and u© = 30 km/s in an orbit inclined at 60° with respect to the galactic plane 


Earth’s revolution around the Sun causes the velocity distribution given in Eq. (3.5), 
and therefore the scattering rate in Eq. (3.1), to modulate in time. Both the time-average 
rate and the modulation amplitude of the rate can be measured. In the SHM, the velocity 
of Earth with respect to the Galaxy is maximum at the end of May or the beginning of 
June. The rate has a maximum at this moment if u m ; n > 200 km/s and it has instead a 
minimum at this moment if u m ; n < 200 km/s. Thus, choosing the modulation phase so 
that the modulation amplitude is positive for u m i n > 200 km/s, the amplitude is negative 
for u m i n < 200 km/s. Only u m i n > 200 km/s are compatible with the phase of the DAMA 
modulation data such that the rate has the maximum on the 2nd of June [7j. The value of 
i^min = 200 km/s is shown as a vertical line in Fig. [I] 

The minimum speed the DM particle must have in the rest frame of the target nuclide 
in order to impart a nuclear recoil energy Er for fiT\d\/m 2 <C 1 is 


^min 


1 

\/2?«tEr 


m T E R 

-bo , 

Ht 


(3.6) 


where jit is the DM-nucleus reduced mass. The mass splitting 6 = m! — m can be either 
positive for endothermic scattering negative for exothermic scattering ED 123 S3 , or 
zero for elastic scattering. By inverting this equation one obtains the minimum and maximum 
recoil energies -Er(u) that are kinematically allowed for a fixed DM speed v, 




(3.7) 


The left panel in Fig. |Tj shows the two E^(v) branches for the Na component of DAMA, 
for 6 = —30 and —50 keV, for WIMP masses that correspond to the best fit regions we will 
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Figure 1. Recoil energy range for exothermic scattering off Na (left) and endothermic scattering 
off I (right), as a function of the WIMP speed v, for the indicated values of WIMP mass m and mass 
splitting 1 5. The two red horizontal lines enclose the 2.0-3.5 keVee energy interval, in which most of 
the DAMA signal is observed. Here Q N a = 0.40 and Q i = 0.09. The DAMA events can only be 
between the two vertical lines at i> m i n = 200 km/s and v m i n = i> max = 765 krn/s (see the text). 


present later. The right panel shows the E^(v) branches for the I component of DAMA, for 
6 = 50 and 100 keV. 

The minimum possible value of v for the interaction to be kinematically allowed is 
vj = \J‘l5 /ht for endothermic scattering, and vj = 0 for exothermic. This speed value 
corresponds to the point of intersection of the two E ^ branches, which occurs at E§ = 
E^(uJ) = E^(vJ) = n T \S\/m T . The maximum value of the DM speed allowed for a given 
halo model is the sum of the escape speed v esc and the modulus of Earth’s velocity in the 
galactic rest frame, which is equal on average to the Sun’s velocity Fq; therefore, for our 
choice of parameter values, w max = 765 km/s. This is indicated with a vertical line in Fig. [lj 
As a result, endothermic scattering with a target nucleus will be kinematically forbidden for 
6 larger than ~ 3.3 keV(/rr/GeV). In the ranges for m and <5 corresponding to the best fit 
regions presented later, endothermic scattering with Na in DAMA is kinematically forbidden. 
On the other hand, scattering off I in DAMA is allowed for endothermic scattering, as can 
be seen in the right panel of Fig. [I] for <5 = 50 and 100 keV. 

In the following sections we examine the compatibility of the WIMP interpretation of 
the DAMA annual modulation signal with various null results for the AV and PS models 
described above. We consider both elastic, endothermic, and exothermic scattering. 

4 Data analysis assuming the SHM 

In this section we describe the data analysis we perform assuming the SHM, which follows 
the procedure already presented in Refs. [241 [48]. 

For the DAMA annual modulation signal, we take the data plotted in Fig. 8 of Ref. [7]. 
We determine the DAMA favored regions in the DM parameter space by performing a Max¬ 
imum Likelihood analysis, assuming the data are Gaussian distributed. Due to the uncer- 
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tainties residing in the quenching factors of Na and I, which play an important role in the 
analysis, we choose two values for each target, namely QNa = 0.40 and 0.30 for Na, and 
Qi = 0.09 and 0.06 [J9| for I (see e.g. Ref. [50] and references therein). In the analysis we 
adopt the combinations Qwa = 0.30 with Qi = 0.06, and Qnb = 0.40 with Q i = 0.09. 

We also compute an upper limit on the WIMP cross section using the total rate measured 
by DAMA, employing the data points plotted in Fig. 1 of Ref. m- We restrict our analysis 
to energies above the experimental threshold of 2 keVee. Given the very large number of 
observed events in each bin, and the resulting small statistical fluctuations, we compute 
an upper bound on the cross section by requiring that the predicted rate does not exceed 
the observed rate in any energy bin. This limit is particularly important for exothermic 
scattering, which reduces the modulation amplitude with respect to the average rate. 

To compute the LUX bound, following Ref. [48] . we apply the Maximum Gap method 
|52] to the variable Si in the range 2-30 photoelectrons. We choose several numbers of 
observed events, i.e. 0, 1, 3, 5 and 24 as described in Ref. [38]. In our SHM model the 
maximum WIMP speed is u max = 765 km/s, thus the maximum recoil energy for a WIMP 
lighter than 11.5 GeV in an elastic scattering with Xe is ~ 12 keV. Using the approximated 
recoil energy contours in Fig. 4 of Ref. [T], and dropping all observed events in and above the 
electron-recoil band (plotted at 1.28 <t) in the same figure, only five observed events remain 
below ~ 12 keV. This means that for elastic scattering, choosing 5 events provides a safe 
upper limit for WIMP masses m < 11.5 GeV. However, if m > 11.5 GeV, only using all the 
24 events lying outside the electron-recoil band provides a reliable upper limit. For inelastic 
scattering with 5 = —50, —30, 50 and 100 keV, the maximum WIMP masses for which 
using the 5 events bound is reliable are 8 . 2 , 9.3, 19.2 and 56.2 GeV, respectively. Since our 
procedure does not depend on the WIMP distribution in the 5i~log 10 (S 2 /S 1 ) plane |T], our 
Maximum Gap upper limits are conservative and safe to be applied to any WIMP-nucleus 
interactions. 

For the SuperCDMS bounds, we use the data set collected by the seven Ge detectors 
between October 2012 and June 2013, corresponding to an effective exposure of 577 kg-days. 
We use the Maximum Gap method with the eleven observed events listed in Table 1 of 
Ref. [2], which passed the selection criteria introduced by the collaboration to discriminate 
signal from background events. We do not incorporate the uncertainty in recoil energy in 
our analysis. For the detector acceptance we take the red curve in Fig. 1 of Ref. [2]. 

To compute the SIMPLE limits, we consider only the Stage 2 [4J, a C 2 CIF 5 detector 
with an exposure of 6.71 kg-day, with one observed event above 8 keV compatible with the 
expected background of 2.2 events. We use the Feldman-Cousins method [53] to place a 90% 
CL upper limit of 2.39 signal events for 2.2 expected background events and 1 observed event. 

For PICASSO we perform a Maximum Likelihood analysis using the data in Fig. 5 
of Ref. [3j. The target material in PICASSO is C 4 F 10 , but the collaboration only considers 
scattering off F in their analysis [3j; we do the same, noting that the contribution of C for DM 
spin-dependent interactions with protons is anyway negligible. We construct our Gaussian 
likelihood using the expected rate above each one of the eight energy thresholds adopted by 
the collaboration (1.7, 2.9, 4.1, 5.8, 6.9, 16.3, 38.8, and 54.8 keV), and the measured rate 
with its uncertainty, which are already background subtracted. 

For KIMS we perform again a Maximum Likelihood analysis using the data points with 
their 68 % CL intervals from Fig. 4 of Ref. [ 6 ], assuming Gaussian distributed data. Because 
Cs and I have similar atomic masses, their quenching factors are not measured separately 
|54j . As for DAMA, we perform our analysis of the KIMS data adopting two values for 
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Q i = Qcs in Csl: 0.05 and 0.10 (see Fig. 2 of [53], and Fig. 5 of [55] and references therein). 

5 Results assuming the SHM 


The plots in Figs. [2]j4j [6]-[8j and 10 13 show 90% CL upper bounds and 68% CL (inner 
and darker shaded region), 90% CL (outer and lighter shaded region), 3cr (solid contour) 
and 5 <t (dashed contour) allowed regions in the m-a p plane. The green shaded regions and 
green closed contours labeled ‘DAMAi’ are the allowed regions compatible with the DAMA 
annual modulation, for quenching factors QNa = 0.40 and Qi = 0.09 in dark green, and 
Q Na = 0.30 and Qi = 0.06 in light green. The lower the quenching factor, the higher is the 
DM mass needed to fit the data. The low and high WIMP mass regions correspond to the 
interpretation of the DAMA data as the WIMP scattering mostly off Na and I in the detector, 
respectively. The upper limit due to the DAMA total rate (black, and labeled ‘DAMAo Na’) 
is shown for scattering off Na assuming QNa = 0.40. 90% CL upper limits from LUX data 
are shown as various magenta curves. As in Ref. [38] the different dashing styles of the lines 
indicate different selections of candidate events used in the Maximum Gap analysis: dotted 
(0 events), double-dot-dashed (1 event), dot-dashed (3 events), dashed (5 events) and solid 
(24 events) curves. Two purple lines show the 90% CL upper limits from KIMS data with 
quenching factors Qi = Qq s = 0.10 (solid) and 0.05 (dashed). 90% CL upper limit from 
SIMPLE (brown), PICASSO (cyan), and SuperCDMS (dark yellow) are also drawn. 



Figure 2. 90% CL bounds and 68% CL, 90% CL, 3er, and 5a allowed regions in the WIMP-proton 
reference cross section a p vs WIMP mass plane, assuming the SHM, for elastic proton-only contact 
AV interactions. The unmodulated DAMA rate limit (black) corresponds to a Na quenching factor 
of 0.40. Different line styles for the LUX bound correspond, from most to least constraining, to 0, 1, 
3, 5 and 24 observed events (see text). The KIMS bound is shown for both Qi = Qc s = 0.10 (solid 
line) and 0.05 (dashed line). 

5.1 Elastic contact interactions 

Fig. | shows our results for elastic proton-only contact AV interactions. The most stringent 
bounds come from two bubble chamber experiments, SIMPLE and PICASSO. Both of these 
limits exclude all the regions favored by the DAMA modulation signal. 
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Figure 3. Same as Fig. [2] but for flavor-universal (left) and proton-only (right) PS interactions. 

Fig. | is the same as Fig. [2j but for PS interactions with flavor-universal coupling 
a n /a p = —1/16.4 (left panel), and with proton-only coupling a n = 0 (right panel). As 
expected, the only limits that change from one case to the other are those of LUX and 
SuperCDMS, due to their enhanced sensitivity to DM-neutron couplings. The DAMA regions 
for WIMP scattering off Na are entirely excluded by SIMPLE and PICASSO, and the regions 
for scattering off I are excluded by KIMS when assuming similar values for the I quenching 
factor in both experiments. This result is different from what was found in Ref. [28], where 
some portion of the Na and I DAMA regions are compatible with all null experiments for 
the PS flavor-universal coupling. 

Ref. [28] uses Bayesian statistics to infer 99% credible level exclusion limits, and 90% 
and 99% credible regions for DAMA, marginalizing over the SHM parameters using Gaussian 
priors (taking central values for the velocities vq = 230 km/s and v = 544 km/s, and for 
the local WIMP density p = 0.3 GeV/cm 3 , with standard deviations Avq = 24.4 km/s, 
Av esc = 39 km/s and A p = 0.13 GeV/cm 3 ). As a result, regions and limits at a specific 
point in parameter space do not necessarily correspond to a fixed set of values for the SHM 
parameters. In our analysis instead we assumed the same set of SHM parameter values across 
all experimental results. We found that the regions and limits move approximately in the 
same manner in the parameter space as we vary the DM velocities, and the DAMA regions 
fail to escape the upper bounds at the 90% CL. This can be seen in Fig. [4] (left panel), which 
shows the results for elastic PS interactions with flavor-universal coupling where both v esc 
and vo are taken 3a below their central values in Ref. [28] (however, we keep u ref = 220 km/s 
to plot a p ). Our choice for the SHM velocities roughly matches the low mass SIMPLE limit 
in Fig. 1 of Ref. [28]. The right panel of Fig. [4]shows also the 99% CL upper bounds (dotted 
lines) for the same set of parameters. In this case, the high mass DAMA region corresponding 
to a quenching factor of 0.09 escapes the KIMS upper limit for quenching factor 0.05. The 
Na component of the DAMA region is still rejected by PICASSO at the 99% CL. In Fig. [5] 
we present our results from the left panels of Figs. [3] and [4] in the same plane as Fig. 1 of 
Ref. |28j, namely in the m-A^ plane, where A^ = m^/y/gjyug. The allowed DAMA regions 
shown in Ref. [28] are much larger than the regions we found. We believe that this is due 
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Figure 4. (left) Same as Fig. [3] but for vq and v esc taken to be 3tr lower than the central values 
assumed in Ref. [2%] . (right) Same as the left panel, but showing in addition 99% CL upper bounds 
(dotted lines) (only the LUX upper bound for 24 events is presented here). 



Figure 5. Same as the left panels in Figs. [3]and|4j but plotted in the = m^/y/gnMg vs WIMP 
mass m plane as in Fig. 1 of Ref. [25] , 


to their marginalization over the SHM parameters and experimental parameters including 
quenching factors. 

Ref. [3TJ found for PS interactions that the LUX bound excludes the DAMA region (see 
Fig. 9 in Ref. [31] . where the I region in DAMA is completely excluded by LUX in both 
the contact and long-range limits). This is in disagreement with our conclusions, possibly 
because of the different analysis of the LUX data. 
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Figure 6. Same as Fig. [2] but for proton-only elastic AV (left) and PS (right) interactions via a 
massless mediator. 


5.2 Elastic long-range interactions 

Fig.[6]shows the regions and limits for elastic AV (left panel) and PS (right panel) interactions 
via a massless mediator. The results are shown in the reference cross section vs DM 
mass m plane, where = a p (m$ = m^ ef ) with = 1 GeV. Note that the results for 
the contact AV and long-range PS interactions are very similar up to a shift in the vertical 
direction (compare Fig. [2] and the right panel of Fig. [6]). This is expected from the Er 
dependence of the differential cross sections given in Eqs. (2.5) and (2.17): disregarding the 


form factors, the differential cross section for the long-range PS and contact AV interactions 
is proportional to E^, for contact PS it is proportional to E^, and for long-range AV it is 
proportional to E p 2 . As it can be seen in Fig. 6j considering long-range elastic interactions 
does not help to bring compatibility between the DAMA regions and the upper limits from 
the experiments with null results. 

5.3 Exothermic contact interactions 

Figs. [7] and [8] show the results for exothermic inelastic proton-only AV and PS interactions, 
respectively, with 6 = —30 keV (left panels) and 5 = —50 keV (right panels). As |<5| increases, 
the DAMA regions move to lower masses, for the following reason. 

The lowest reach in DM mass for a direct detection experiment is obtained when 
i = (see Fig. [lj), where the threshold energy E t ^ is the lowest detectable re¬ 

coil energy. The mass reach can be found by extracting m as a function of 5 and v from 
Eq. ( |3T| ), for E+(v) = E th , 


Er_( v i 


? b (< 5 , v ) = 


E th m T 


vy/2E th m T - E th - 5 ’ 


(5.1) 


and evaluating it at v = u max . The DAMA region will be therefore located at WIMP masses 
higher than fh(5,v ma , x ), taking Na as the target element. For DAMA, E t h = 5 keV for 
scattering off Na with quenching factor 0.40. The lowest reach of DAMA is the lower green 
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Figure 7. Same as Fig. [2] but for exothermic AV interactions with S = —30 keV (left) and S = —50 
keV (right). 
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Figure 8. Same as Fig. [7] but for PS interactions. 


line plotted in Fig. [9j Also shown in the same figure are the mass reaches of PICASSO and 
SIMPLE, for which we used E t h = 1.7 keV and 8 keV, respectively. We only considered 
scattering off F in both experiments. 

An estimated upper limit on the DM mass for the DAMA region comes from requir¬ 
ing that DM particles with speeds below 200 km/s always scatter below threshold and are 
therefore undetectable, because otherwise DAMA should have observed a sign change in the 
modulation amplitude. In other words, scatterings of DM particles slower than 200 km/s 
would yield a different phase for the modulated signal with respect to that measured by 
DAMA, and therefore an acceptable fit requires these scatterings to occur below threshold. 
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Figure 9. Estimated lowest reach in WIMP mass m for SIMPLE and PICASSO, and estimated 
mass range in which the Na component of the DAM A region with quenching factor 0.40 is found, as 
a function of mass splitting S for exothermic scattering. 


Since a fixed nuclear recoil energy can be imparted by heavier DM particles traveling at lower 
velocities, the condition n m i n (^th) > 200 km/s implies an upper limit on the DM mass in 
DAMA, given by rh(5,v = 200 km/s). This upper limit is the higher green line plotted in 
Fig. i The other possible condition to avoid scatterings of DM particles slower than 200 
km/s, i.e. having a large enough E^(v = 200 km/s), would imply a very odd spectrum in 
DAMA, with more events at higher energy instead of the observed spectrum vanishing at 
high energy. 

Since exothermic scattering decreases the value of u m i n for a given recoil energy, the 
modulation amplitude becomes smaller with respect to the time-average rate. For large 
enough |<5| the DAMA modulation signal becomes inconsistent with the DAMA time-average 
rate. For values of 5 lower than about —30 keV (for AV interactions) and —50 keV (for PS 
interactions), the DAMA total rate limit rules out the modulation signal in Na, as indicated 
by the black curve excluding the DAMA Na region in the right panels of Figs. [7]and[8j 

For the values of 5 allowed by the DAMA rate, the limit by PICASSO (and also the 
SIMPLE limit in most instances) rejects the allowed regions. For each value of 6 on the 
horizontal axis of Fig. [9j the DAMA region spans a mass range enclosed within the green 
belt, while the SIMPLE and PICASSO lines indicate the mass value where the limits in the 
m-(jp plane become vertical. From the plot it becomes clear that exothermic scattering brings 
compatibility between SIMPLE and DAMA for large enough |<5|, as suggested by Figs. [T] and 
[8j however the region is rejected by the DAMA average rate measurements. While this is true 
for QNa = 0.40, smaller quenching factors move the DAMA region to larger DM masses, thus 
potentially compromising this compatibility with SIMPLE. In any case, the DAMA region 
does not escape the PICASSO limit. 
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5.4 Exothermic long-range interactions 

Fig. [TO] shows the AV interaction via a massless mediator for 5 = —30 keV. Here as well, all 
the DAM A regions are rejected by the null experiments. We do not plot the results for long- 
range PS interactions as these are qualitatively similar to those for contact AV interactions, 
as commented above. 



Figure 10. Same as Fig. [7] (left) but for a massless mediator. 


5.5 Endothermic contact interactions 


Figs. m and [12] show the result for the proton-only spin dependent endothermic scattering 
with AV and PS interactions, respectively, in the contact limit, with 6 = 50 keV (left panels) 
and 5 = 100 keV (right panels). As 5 increases, scattering off light targets becomes kinemat¬ 
ically forbidden since vj becomes larger than u max . For 6 = 100 keV, the only remaining 
limits are from KIMS and LUX. We can see that the DAMA region for I scattering moves 
towards the left compared to the KIMS upper bound as 5 increases. For PS interactions 
with 5 = 50 keV, it is only the combination of larger quenching factor Qi = 0.09 for I in 
DAMA and smaller quenching factor Q\ = 0.05 in KIMS that allows the DAMA signal to be 
compatible with all upper limits. For 6 = 100 keV, the I region corresponding to Qi = 0.06 
barely escapes the limit with Q 1 = 0.05 from KIMS, and the situation remains tense for 
quenching factors 0.09 (DAMA) and 0.10 (KIMS). Raising 5 further makes it progressively 
more difficult to find a region of the DAMA signal that is kinematically allowed. For AV in¬ 
teractions the DAMA regions are even more severely constrained: only the larger quenching 
factor Qi = 0.09 DAMA region is allowed by the KIMS upper bound with smaller quenching 
factor Qi = 0.05 for 6 = 100 keV. 

These results are largely consistent with those of Ref. |20| , where the framework of 
non-relativistic operators introduced in Ref. [32] was generalized to inelastic scattering and 
a model-independent analysis was performed on a series of effective operators. In Table V of 
Ref. [20] the AV and PS interactions correspond to fermion operators 15 and 4, respectively. 
For those interactions, Ref. [20] quotes best fit parameters for DAMA (corresponding to 
5 = 106 keV for AV and S = 57 keV for PS) that are consistent with the KIMS data only 
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Figure 11. Same as Fig. [2] but for endothermic AV interactions with 8 = 50 keV (left) and <5 = 100 
keV (right). 




Figure 12. Same as Fig. EH but for PS interactions. 


for DAMA quenching factor 0.09 and KIMS quenching factor 0.05. In Ref. [20], however, 
scattering off Cs in KIMS was neglected due to the lack of the form factor in Ref. [32] . 
Since the contribution of Cs to the scattering rate is sizable for interaction with protons, we 
adopt an approximate form factor for Cs as discussed in Sec. 2.3 resulting in stronger KIMS 
bounds. 

At this point it is important to recall the flavor physics bounds we mentioned in the 
introduction. Fig. 9 of Ref. m shows that the quark couplings needed for PS inelastic 
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scattering on I to fit the DAMA data with a one-particle exchange process, 


9 5dm = 


2 

V100 MeV/ 


for 100 MeV 

for <C 100 MeV, 


(5.2) 


where k — 0.2 and k! — 0.1 for S = 50 keV, and k ~ 1.6 and k! ~ 0.7 for <5 = 100 keV, 
is rejected m for any reasonable value of <?dm (unless #dm > 10 5 ). Note that the quark 
coupling used in Ref. m is g/V 2. 


5.6 Endothermic long-range interactions 


For the AV interaction via a massless mediator shown in Fig. 13, only a very small portion 


of the I DAMA region with quenching factor 0.09 escapes the KIMS limit with quenching 
factor 0.05. Therefore, quite different quenching factors are needed for I in DAMA and KIMS 
to have the DAMA region escape the KIMS limit. Although quenching factors of a given 
element in different crystals can have different values in general, large differences may be 
questionable. Again, we do not plot the results for long-range PS interactions since these are 
qualitatively similar to those for contact AV interactions, as commented above. 



Figure 13. 


Same as Fig. 11 (right) but for a massless mediator. 


6 Halo-independent analysis 


So far we have assumed a particular model for the dark halo of our galaxy. It is however 
possible to compare direct detection data without making any assumption about the local 
density or velocity distribution of the dark matter particles US \22l SHI [56H75] (in particular 
we follow the analysis of Refs. da hhi E6 i E7]). This method consists of extracting from the 
data, instead of just the (reference) WIMP-proton cross section a p , the function fj of v m \ n 
defined below in Eq. (6.2) which encloses all the dependence of the rate on the DM velocity 
distribution. Since this function is experiment-independent, data from different experiments 
can be directly compared in the v m \ n -'q plane. 
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In order to perform the halo-independent analysis we have to assume a value for the 
DM mass to, together with all the other interaction parameters such as the mass splitting 
5 and the neutron to proton coupling ratio a n /a p . We study parameter values that seem 
promising in our SHM analysis to make a DM interpretation of the DAMA data compatible 
with all other experiments when relaxing the assumption on the dark halo. Taking into 
account Figs. [3] and [4] we select a WIMP mass below 10 GeV, one close to 30 GeV, and one 
close to 50 GeV, for WIMPs with PS interactions and elastic scattering. For the couplings 
we take a n = 0 rather than a n /a p = —1/16.4; this choice is conservative in the sense that, 
while results from experiments employing targets with negligible spin-dependent interactions 
with neutrons like DAMA, SIMPLE, PICASSO, and KIMS are not affected, the bounds 
from LUX and SuperCDMS are less constraining when the WIMP-neutron coupling is set to 
0 (see e.g. Fig. [3]). We do not consider inelastic exothermic scattering, as the DAMA regions 
are badly excluded in all cases studied in the previous section. For inelastic endothermic 


scattering, looking at Fig. 11 we select m = 40 GeV for 5 = 50 keV and m = 52 GeV for 


5 = 100 keV, for WIMPs with AV interactions with a n = 0. Analogously, from Fig. 12 we 
select to = 38 GeV for <5 = 50 keV and m = 45 GeV for <5 = 100 keV, for WIMPs with 
PS interactions with a n = 0. Notice that some of these choices are similar to the best fit 
parameters of Ref. [20], i.e. to = 54.3 GeV and <5 = 106 keV for AV interactions and to = 40.8 
GeV and 5 = 57 keV for PS interactions. Finally, from Fig. [13] we select to = 80 GeV for 

6 = 100 keV, for long-range AV interactions with a n = 0. 


6.1 More on the direct detection rate 


Since the differential cross section dax/dE R for the AV and PS interactions is proportional 
to 1/v 2 (see section]^, and thus v 2 dax/dE R is independent of v, we can express Eq. (3.2) as 


R[E[,E^t) = J2— rdE R fj( Viain (ER),t)(-^f-) f^ 2 dE 1 e(E\ E R )G T (E R , E') , 

rj-1 r fn r [' J Q \ R / i/ 

(6-1) 

where a p is the (reference) total DM-proton scattering cross section and we defined the 
velocity integral 

rj(v m \ n ,t) = — [ d 3 v _ (6.2) 


Notice that, while fj as a function of v m \ n is independent of the target nuclide, v m \n {E R ) is a 
function of E R that depends on the target. 


By changing the integration variable from E R to v m i n in Eq. (6.1), the event rate can 
be written in the more compact form 


POO 

= 'y ^ / du m i n ?/(e m i n ,f)(R[g/ g/j(u m i n ) , 
rn J vT 


(6.3) 
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with vj = \J25 /ht. The target-specific response function OlJ E , £,j(u m i n ) is given by 
■^[e' 1 z,E' 2 ]( V min) = 
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dry 
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(6.4) 


with E R (v m i n ) given in Eq. (3.7). For elastic scattering (6 = 0), E R = 0 and E R equals 


ER{v min ,m T ) = ‘ 2 ^ L vl l] 
m.T 


(6.5) 


so that IRj^g, E /|(f min ) reduces to 


^fE;, E 2]( t ’ min ) — 


Ct 4 liLv- 
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rriT rriT 


v 2 d ax 
cr p dE R 


(E R (v min ,m T ),v) 


rK 

/ d E' e(E' 

Je[ 


,E R (v m i n ,m T ))G T (E R (v m i n ,m T ),E') . (6.6) 


Due to Earth’s rotation around the Sun, the velocity integral (6.2) is modulated in time 
with a 1 year period: 


iK'Dnin ,t) ^ fj°(v m in) + fj 1 (v min ) COS ( — (t - t 0 ) ) . 

\y r J 

Since all the time dependence of the rate is contained in fj, we also have that 

^'[E' 1 ,E' 2 ]('t) — -R[ £',£'] + R[E[,E' 2 ] cos ~ *o)) J 


(6.7) 


( 6 . 8 ) 


where to is the time when the rate reaches its maximum value, and R? E , E ,, and Ri E , E ,^ are 
the time-average and annual modulation amplitude of the rate, respectively. In general fj can 
be expanded in a Fourier series, and here we assume that higher modes are not important. 
The DAMA collaboration did not find any hints of higher modes in their data [71 [76], thus 
when considering the DAMA data we adopt Eq. (6.7) with the measured phase to = June 
2 nd . All other experiments considered here give an upper bound on a time-averaged signal, 
thus on fj°. Given that the annual modulation amplitude cannot be larger than the average 
rate, |r) 1 (u T nin )| ^ f?°(u m i n ), we can interpret upper bounds on fj° from experiments with null 
results as (conservative) limits on the fj 1 signal measured by DAMA. We do not consider the 
direct CDMS-II bound on fj since LUX and SuperCDMS set more stringent constraints, 
see Refs. [13, 241148] (notice also that SuperCDMS employs the same target material). 

6.2 Data analysis for the halo-independent method 

For LUX, SuperCDMS, and SIMPLE we follow the procedure developed and described in 
Refs. [T3l I4H1 1561 157] . For PICASSO and KIMS we cannot perform a Maximum Gap analysis 
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as done for LUX and SuperCDMS because the data are binned. We therefore produce a limit 
on fj° at each U min value in the following way. We compute the rate (6.3) adopting a step 
function, 77°(ti min ') = fj*9( — u m i n ), because it is the function that allows to draw the most 
conservative bound on the value fj* taken by ? 7 0 (u m i n ) at a specific u m i n value v U n [521 RTTI . 
For each value of iW n we compare this predicted rate in each single energy bin with the 90% 
CL limit on the rate. Imposing the computed rate not to surpass the limit in any of the bins 
thus fixes the maximum allowed fj* at ujL n . For KIMS we use the black limit lines in Fig. 4 of 
Ref. [6j, while for PICASSO we translate the upper end of the error bars in Fig. 5 of Ref. [3] 
into 90% CL upper limits assuming the data are Gaussian distributed and the uncertainty is 
given at the la level. 

The halo-independent analysis of the DAMA annual modulation data presented in 
Refs. 113 , i S51 [56l I5T] and in Ref. [22] is only applicable when WIMPs can scatter off only 
one of the target elements, either Na or I. This happens, for instance, if the DM is so light 
that elastic scattering off I occurs always below threshold, assuming a reasonable maximum 
speed with respect to Earth, u max , for WIMPs in the galaxy. It also happens for inelastic 
endothermic scattering with u max < = \j2b //7Na, which makes WIMP scattering off Na 

kinematically forbidden. Therefore, we can straightforwardly apply the analysis of DAMA 
data presented in Refs. [13, 22 HU -16 22 to the light WIMP with m < 10 GeV scattering 
elastically and to all considered cases of endothermic scattering. The WIMPs with mass 
close to 30 GeV and 50 GeV scattering elastically need a special treatment, described in the 
following. 

When both elements, Na and I, are involved in the scattering, the same value of the 
detected energy E' is mapped onto different u m i n values because of the different target masses 
and quenching factors. In order to extract the value of fj in a u m i n interval from Eq. (6.3), 
we adapt to the DAMA data the procedure that was developed for CRESST-II in Appendix 
A.2 of Ref. [M] an d in Ref. [23] • We start by choosing the highest energy bin at high enough 
energy so that the heaviest element (I) does not contribute to the scattering rate in it, 
because the necessary u m i n exceeds the maximum possible speed in the halo. Starting from 
the highest energy bin we work our way toward the lowest energy bin as depicted in Fig. [Lf| 
We compute the u m i n range corresponding to the highest E' bin for Na using the relation 
(E') = QTE-R(v m i n , 777 ,r) for T = Na (with En(v m i n ,mT) given in Eq. (6.5)). From this 7; m i n 
range we derive the second highest energy bin using again the same relation for T = I. The 
procedure can then be repeated starting from this new bin, until the lightest energy bin above 
the experimental threshold is built. Since we only wish to consider bins where a significant 
signal is observed, we require the lowest energy bin to lie within the 2.0-6.0 keVee interval, 
where the modulation amplitude measured by DAMA is significantly different from zero. 

Notice that the DAMA modulation data span 2.0 to 20.0 keVee with an original bin size 
of 0.5 keVee. In rebinning the DAMA data we want to merge the original bins, and not split 
existing bins. Therefore, each new bin needs to have boundaries that are multiples of 0.5 
keVee so as to encompass an entire number of original DAMA bins. Notice also that, if E' is 
the boundary of one of the chosen bins, the corresponding boundary of the next higher energy 
bin is rE' with r = QN a E R (v min ,rn Na ,)/QiE R {v min ,rni) = QNa^Na m i/ Qn4 m Na.- Therefore, 
it is necessary for r to be an integer (or a half-integer) number. For quenching factors that 
are constant in energy, this can be achieved for particular DM mass values. With Qn a = 0.30 
and Qi = 0.09, we choose our DM particle masses to be 30.14 GeV and 47.35 GeV for which 
r = 5.0 and 3.5, respectively. 

Choosing u max = 800 km/s, for m = 30.14 GeV the bin at highest E' must be completely 
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DAMA data 





Figure 14. Binning scheme used for m = 30.14 GeV and elastic scattering. The red and blue 
arrows show the correspondence between average detected energies ( E '), shown on the top axis, and 
minimum speeds c m i n , shown on the bottom axis, for Na and I targets, respectively. The dashed 
arrows show that WIMPs scattering off I could produce events in the highest energy bin only with 
c rn in > c max = 800 km/s. The DAMA data are overlaid in green. 


above QtEk mi) — 6.3 keVee for Qi = 0.09. Using r = 5 we take two bins, [2.0,4.0] 
keVee and [10.0, 20.0] keVee. We choose bins well separated in energy and as large as possible 
to avoid overlapping and to minimize the effect of the tails of the corresponding response 
functions. The binning scheme for this candidate is shown in Fig. 14 Our choice of m = 47.35 
GeV comes from a halo model in Fig. [4] with a low u max value, close to u max = 600 km/s. 
Assuming this u max value, scattering of I is kinematically forbidden for E' energies above 
Q.jE-r (u m a y , tut) ~ 6.97 keVee for Q\ = 0.09. Using r = 3 we choose the two bins [3.0,6.0] 
keVee and [10.5, 21.0] keVee. Since the highest energy bin surpasses the 2.0-20.0 keVee energy 
range where the DAMA modulation data is available, for the two additional 0.5 keVee bins 
in the 20.0 to 21.0 keVee range we assume the same average and mean square error as for 
the nineteen 0.5 keVee bins in the 10.5-20.0 keVee range. 

Once the energy bins to be used in the analysis are established, we extract information 
on the modulated component of the velocity integral, fj 1 , in the following way. For two bins 


[E[, E' 2 \ and [E^, E'^] (the extension to a larger number of bins is trivial), Eq. (6.3) reads 


R 


re 

[E[,E^] = J dVminfj (Unin) (^min) + (^min) 


= A 


Na 


1, Na 


[E'i,E^\V [E [,E 


+ A [E' 1 ,E' 2 ] f i[E[,E! l 


(6.9) 


(Unin) 
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R 


[E' 3 ,Ei] 


du m i n f) (u m i n ) ,E’^y min) T -^-[E 1 ,E' 4 ] (^min) 


~ A 


Na 


1, Na 


[e' 3 ,e' 4 ]V [E ' 3 ,e 


( 6 . 10 ) 


where scattering off I does not contribute to the rate in the highest energy bin by construction. 
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Figure 15. Normalized response functions 31 T (u m in) /-A T for WIMPs with PS interactions and elastic 
scattering, for DM mass m = 30.14 GeV (left) and m = 47.35 GeV (right). 1R|^, £ v](timin)A4[£' is 
shown in red and 1R| S , ^(^mi n )/A^ E , E ,^ in blue. 


[E[,E£ 


Here we defined the target-specific average of fj^ E , E ,j as 

/ 0 °° dn min fj(v min )We;]( ^min) 


~1 ,T _ 


a t 


with 


/»oo 

= J dUmin lR[ B / iE /](^ m i n ) , 


( 6 . 11 ) 


( 6 . 12 ) 


and analogous definitions for [Eg, £ 4 ]. Since the two energy bins are chosen so that scattering 
off I in [E[, E' 2 \ probes the same u m ; n range as scattering off Na in [£ 3 , £ 4 ], we expect that 


~l,Na ~1,I 


(6.13) 


and thus that E ,](u min )/.Aj|| E ,^ ~ 3?j E , E ^{v min )/A 1 ^ E ,^. Fig. |l5| shows the comparison 
between these quantities for WIMPs with PS interactions and elastic scattering, and gives 
an idea of the extent to which the assumption in Eq. (6.13) is correct. In Fig. 15 we can 
see that the v m i n range associated with a certain detected energy bin through the average 
relation (£') = QtE-r (v m\n , mr) determines only approximately the u m i n range in which 
the corresponding response functions ^^(umin) are significantly different from zero. The 
difference between the response functions of Na and I shown in Fig. [15] is due to the width of 
the energy resolution function, which depends on energy, and the u m i n (or alternatively the 
£r) dependence of the scattering cross section. 


Finally, using Eqs. (6.9), (6.10) and (6.13) we get the system of equations 


pl 
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_ /jNa Na 

^[E' x ,E'^[E^E' 2 


+ “df rpt 


~ 1 ,I 
\E[,E^\E[,E ' 2 


(6.14) 


R 


/iNa ~1, Na 
[^ 3 .^ 4 ] - A [E' 3 ,Erie's,E' 4 ] 


~M Na 


~ 1 ,I 

[ E'^E'/llE'^ 


- 24 - 



















Figure 16. 90% CL upper bounds on 77°(u m i n ) from LUX, SuperCDMS, SIMPLE, PICASSO, and 
KIMS, and measurements of r? 1 (u m ; n ) with 68% CL vertical error bars for DAMA with Q^ a = 0.30, 
Qi = 0.09, for a WIMP with mass m = 7 GeV (top), to = 30.14 GeV (bottom left) and 47.35 GeV 
(bottom right), contact PS interactions and elastic scattering. Different line styles for the LUX bound 
correspond, from most constraining to least constraining, to 0, 1, 3, 5 and 24 observed events (see 
Section|4|. The KIMS bound is shown for both Q\ = Qc s = 0.10 (solid line) and 0.05 (dashed line). 
The thin DAMA crosses show the absolute value of fj 1 when this is negative. 


which can be solved for an d v\e' E'y We com P u te the 68% CL error on these quantities 

by propagating the la uncertainty of the DAMA data. 


6.3 Results of the halo-independent data comparison 

The plots in Figs. 16-19 show 90% CL upper bounds on fj°(v m j n ) from LUX, SuperCDMS, 
SIMPLE, PICASSO, and KIMS with Qi = Qq s = 0.10 (solid purple line) and 0.05 (dashed 
purple line). The DAMA measurements of the annual modulation amplitude f) 1 (u m i n ) are 
shown as crosses, where the vertical bars show the 68% CL uncertainty and are located at 
the position of the maximum of the relevant response functions 3? r (u m i n ). The horizontal 
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bar of each cross indicates the w m i n interval where 90% of the integral of lk T (w m i n ) about the 
peak is included. When both Na and I contribute to a cross we choose the response function 
extended over the larger w m i n interval, which is J^fu min ) (see Fig. 15). We assume the most 
commonly adopted values of the quenching factors, = 0-30 and Qi = 0.09. 

For WIMPs with PS interactions and elastic scattering (see Fig. 16), we selected three 
masses from the DAM A regions shown in Figs. [3] and [4j We show results for a n = 0, as 
explained above. Only scattering off Na is kinematically accessible in DAMA for m = 7 GeV 
(top panel of Fig. 16), since E' > 2 keVee would require v m \ n > 1644 km/s for I recoils. In this 
case the limits of PICASSO and SIMPLE cut across the DAMA points (each corresponding 
to the DAMA bins of width 0.5 keVee from 2.0 to 6.5 keVee), except for the highest energy 
bin. This shows incompatibility between the DAMA and the PICASSO and SIMPLE data, 
unless the modulation amplitude |f) 1 (u m i n )| can be as large as the time-average fj°(v min ), 
which is not possible in the whole u m i n range. 

For the bottom panels of Fig. 16, where m = 30.14 GeV and 47.35 GeV, both I and Na 
contribute (see above). Although these masses were chosen within the DAMA regions in our 
SHM analysis, the fj 1 crosses shown in the two panels correspond to a halo that significantly 
differs from the SHM. We obtained negative fj 1 values in the lowest v m i n ranges in both 
cases, both of them above u m i n = 200 km/s, while in the SHM r) l is negative only below the 
Dnin = 200 km/s value. The absolute value of the crosses with negative fj 1 values are shown 
in Fig. 


16 


' with thinner lines than those with positive fj 1 values. In both cases the combination 
of KIMS, PICASSO, SIMPLE and LUX bounds reject the DAMA crosses, showing strong 
incompatibility between the DAMA result and the just mentioned limits. 

Our results are similar to those of Ref. m, although the halo-independent analysis 
done in this reference is different. In Ref. m the “minimal” f) 1 (n m i n ) compatible with the 
lcr DAMA error bars is identified with the piecewise continuous function touching the lower 
end of these error bars. Only the four 0.5 keVee experimental bins in the 2.0 to 4.0 keVee 
energy interval are considered. ? 7 1 (n m i n ) is set to zero outside the 'u m j n range corresponding 
to this E' range. This amounts to a choice of very low (m-dependent) u max . With this 
choice, the contribution of I can be neglected for m < 60 GeV. Besides, ? 7 0 (u m i n ) is assumed 
to be equal to the “minimal” ^ 1 (n m i n ) in almost all the considered u m i n bins. This implies 
the dark halo leading to this velocity integral is quite different from the SHM, as we also 
find. However, while our results seem similar to those of Ref. m, we partially draw different 
conclusions. The authors of Ref. m conservatively conclude that their choice of fj 1 and fj° 
for DAMA can be compatible with the limits from other experiments for PS interactions 
with m = 7 and m = 30 GeV (while being rejected for AV interactions). This is because 
they allow r) 1 (n m i n ) to be equal to ?j°(u m i n ) in almost the whole 2.0-4.0 keVee energy range. 
Assuming instead that |fj 1 (v m i n )| is much smaller than ij°(u m i n ), as it happens in most halo 
models, the DM interpretation of the DAMA data is in strong tension with the limits. 

Fig. 17 shows the result of our halo-independent analysis for WIMPs with AV interac¬ 
tions and inelastic endothermic scattering. The left panel is for m = 40 GeV and 5 = 50 
keV, for which i?? a = 802.7 km/s, and the right panel is for m = 52 GeV and 5 = 100 
keV, for which = 1135.2 km/s. Therefore, only scattering off I is kinematically allowed 
for Umax = 800 km/s. The DAMA data seems in disagreement with the KIMS bound for 
Q i = Qcs = 0.10, but not if Q\ = Qcs = 0.05 in KIMS (however, this last value is much 
smaller than the Q\ = 0.09 taken for DAMA). 

Fig. [18] shows the result of our halo-independent analysis for WIMPs with PS inter¬ 
actions and inelastic endothermic scattering, for m = 38 GeV and 5 = 50 keV in the left 
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v mi „ [km/s] v min [km/s] 


Figure 17. Same as Fig. [l6|but for WIMPs with inelastic endothermic contact AV interactions and 
to = 58 GeV, S = 50 keV (left), and to = 52 GeV, S = 100 keV (right). 




v ml „ [km/s] 


v mln [km/s] 


Figure 18. Same as Fig. [IT] but for WIMPs with inelastic endothermic contact PS interactions and 
to = 38 GeV, S = 50 keV (left), and to = 45 GeV, S = 100 keV (right). 


panel, resulting in = 810.1 km/s, and m = 45 GeV and 6 = 100 keV in the right panel, 
resulting in = 1113.0 km/s. WIMPs scattering off Na are again kinematically forbidden 
for u max = 800 km/s. The DAMA data are in disagreement with the KIMS bound, unless 
Qi = Qc s = 0.05 and |i/ 1 ( ?; min)| - fj°(v min ). 

Fig. [19] shows our results for a WIMP with AV interactions and long-range inelastic 
endothermic scattering. Here m = 80 GeV and 5 = 100 keV, thus r// a = 1031.5 km/s. The 
DAMA data are excluded by our most stringent KIMS (Qi = Qcs = 0.10) and LUX bounds, 
but they remain in strong tension with the null results even when the less stringent KIMS 
and LUX bounds are considered. 
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Figure 19. Same as Fig. [IT] but for WIMPs with inelastic endothermic AV long-range interactions 
and to = 80 GeV, S = 100 keV. 


7 Conclusions 


We investigated the possibility of interpreting the annual modulation signal observed in 
the DAMA experiment as due to WIMPs with spin-dependent coupling mostly to protons. 
We considered both an axial-vector (AV) interaction, which is what is usually referred to 
as ‘spin-dependent interaction’, and a pseudo-scalar (PS) interaction, proposed in Ref. |28] 
to reconcile DAMA with the null experiments. We also extended our analysis to inelastic 
scattering, and considered both contact and long-range interactions. Due to the similar Er 
dependence of the differential cross sections, we find for the long-range PS interaction the 
same results as for contact AV interactions, up to a shift in a p . We analyzed the data both 
assuming the Standard Halo Model (SHM) and in a halo-independent manner. 

Spin-dependent WIMP couplings mostly to protons effectively weaken the bounds from 
experiments using Xe and Ge as target elements, whose spin is due mostly to neutrons. 
However, the bounds from experiments with F and I targets, such as PICASSO, SIMPLE 
and KIMS, remain relevant since their spin is due mostly to protons. 

Assuming the SHM, for elastic scattering (see Figs. [2] to [6]) we found that, in all the 
cases investigated here, the DAMA regions for Na are entirely excluded by SIMPLE and 
PICASSO, while the regions for I are excluded by KIMS. 

For exothermic scattering (see Figs. 0® and [To| , the DAMA regions move progres¬ 
sively to smaller WIMP masses with respect to the upper limits as |<5| increases, because the 
modulation phase observed by DAMA forces u m i n > 200 km/s, and this is possible only for 
progressively lighter WIMPs (see Fig. [9]). Thus, exothermic scattering brings compatibility 
between the DAMA region for Na and the upper bounds from SIMPLE. However, it does not 
suppress the PICASSO limit which continues to rule out the DAMA region. Furthermore, 
exothermic scattering reduces the modulation amplitude with respect to the time-average 
rate. Thus, the upper limit derived from the DAMA average rate measurement rejects the 
interpretation of the signal as due to scattering off Na for values of 6 < —30 keV. The DAMA 
region for scattering off I is excluded by the SIMPLE and KIMS upper bounds. 
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For endothermic scattering (see Figs. 11 to 13), only KIMS provides relevant bounds. 


Scattering in all detectors besides KIMS and LUX becomes kinematically forbidden for large 
enough S. We showed results for <5 = 50 and 100 keV, because as 5 increases further, scattering 
off I becomes kinematically forbidden as well. For 6 = 50 keV, only assuming a larger 
quenching factor Q\ = 0.09 for I in DAMA, and a smaller quenching factor Q\ = Qq s = 0.05 
in KIMS, the allowed DAMA region is compatible with all present limits for PS couplings. 
However, the possibility that the same nuclide has such different quenching factors in different 
crystals may be questionable. The same holds for contact and long-range AV and long-range 
PS interactions for 5 = 100 keV. For contact PS interactions and 5 = 100 keV, a small sleeve 
of the 90% CL DAMA region for scattering off I escapes the 90% CL KIMS limit with similar 
Qi for both experiments. These results are largely consistent with the results of Ref. |20j . 
However, for PS interactions the DAMA regions are rejected by flavor physics bounds on 
the PS coupling to quarks [31] (unless ^dm can be very large (?dm > 10 5 ). In our analysis 
we assumed that the scattering process can be approximated by one-particle exchange. The 
inclusion of multi-particle exchange processes may change the form of the WIMP-nucleus 
scattering cross section, and therefore all bounds shall be reconsidered. 

For WIMP mass values within the DAMA regions derived assuming the SHM and close 
to the upper limits rejecting them, we performed also a halo-independent analysis. This is 
presented in Fig. 16 for elastic scattering with contact PS interactions and in Figs. 17-19 for 
inelastic endothermic scattering with contact AV, contact PS, and long-range AV interactions. 
We again find strong tension between the DAMA data and upper bounds, except for contact 
AV and PS interactions with inelastic endothermic scattering if Q\ in KIMS is much smaller 
( Q\ = 0.05) than the Q\ = 0.09 assumed for DAMA, although this choice of different Q\ 
values for both experiments may be questionable. 
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